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Abstract 

The dynamics of the discrete Gaussian model for the surface of a crystal 
deposited on a disordered substrate is investigated by Monte Carlo simula- 
tions. The mobility of the growing surface was studied as a function of a small 
driving force F and temperature T. A continuous transition is found from 
high-temperature phase characterized by linear response to a low-temperature 
phase with nonlinear, temperature dependent response. In the simulated 
regime of driving force the numerical results are in general agreement with 

recent dynamic renormalization group predictions. 
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There has been considerable progress in recent investigations of crystalline surface growth 
It is known [||,[| that, due to the discreteness and fluctuations (thermal fluctuations 
and fluctuations in the growth process itself), a crystalline surface undergoes a phase tran- 
sition between a high-temperature rough phase and a low-temperature smooth phase. The 
presence of the quenched disorder in the crystal (either in the substrate [[]] or in the bulk 
0) changes both the critical temperature and the character of the low-temperature phase. 
Below critical temperature, T c , the system develops a glassy phase characterized by the 
existence of many metastable states to which surface configurations are pinned by disorder. 
The surface itself remains rough but with very different static and dynamic properties ||[7] 
compared to the roughness above T c . 

In the present work the dynamic properties of the roughening transition are examined 
by numerical simulations. The numerical studies of the static properties were reported 
elsewhere ||. The simulated system is based on the Hamiltonian of the discrete Gaussian 
model which was very successful |||| in describing the surface on a flat substrate: 

H = £ £ ~ h 3 f. (1) 

The sum runs over nearest-neighbor pairs, k is the surface tension, and hi is the height of 
the surface above the point % on the two-dimensional basal lattice. In the case of a flat 
surface hi takes integer values in the units of lattice spacing a in the direction perpendicular 
to the surface. To simulate the disordered substrate a random quenched height di chosen 
uniformly (and independently for each site) in the interval (— a/2, +a/2] was first assigned 
to each site. The height hi then takes the values hi = di + n^a where rii is any positive or 
negative integer. 

In the continuum limit, hi — > 4>{x), di — > d(x), and near the critical point, Eq. ([]]) maps 
to the Hamiltonian of the random phase sine-Gordon model (RSGM) [|J: 

H = [V0(f)f - g cos (2tt [0(f) - d(x)\ /a) j . (2) 

The periodic cosine term comes from the lattice discreteness and is crucial for the exis- 
tence of a phase transition. The constant g might be considered as the strength of the 
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periodic pinning potential. The Hamiltonian (0) also describes other disordered systems: 
two-dimensional vortex glass with a parallel magnetic field |9|], charge density waves, and it 
is also equivalent to the vortex-free XY model with random field. 

If the disorder is absent, the predicted value [|IJ for the critical temperature of the discrete 
Gaussian model ([]]) is of the order Tr ~ 1.45, which is close to the results from computer 
simulations [fLOfl . The behavior of the roughening transition in the presence of an applied 
force F was studied in detail by Nozieres and Gallet (NG) ||. They found a broadening 
of the transition due to nonequilibrium conditions with crossover occurring below Tr. In 
the limit F — > 0, the interface mobility fi = v/F has a sharp jump at Tr, from a finite 
value for T > Tr to zero for T < Tr. If surface is driven by a small, but finite, driving 
force it remains rough (with temperature- and force-dependent mobility) even below Tr. 
The NG theory describing the system in the temperature region close to Tr ("interrupted 
renormalization" scheme) or well below Tr (homogeneous nucleation) is in quantitative 
agreement with experiments on (0001) interfaces of hep 4 He crystals near the roughening 
transition [IT]. 



The effects of the disordered substrate on the dynamic properties in the vicinity of the 
phase transition were studied by the dynamic renormalization group (RG) methods p.|T2 
based on the Langevin dynamics with a small driving force F: 

Here the Gaussian fluctuating noise £ satisfies t)£(x', 0)) = 2T5(x — x')5(t) and 
the exponentiated random phase d(x) in Eq. ($) also obeys Gaussian statistics: 
(exp{t27cd(x)/a} exp{— t2nd(x')/a}) = a 2 5(x — x'). The theory predicts that above the 
critical temperature, T c = k/tt, the system responses linearly to the applied force, i.e., 
the mobility fi is finite and independent of F, while below T c the response is nonlinear, 
characterized by the temperature-dependent exponent r\: 

(T/T c -1)< for T>T C , 
V- ~ { (4) 
(1 - T/TcYF* for T < T c , 



where £ « 1.78 is a universal constant and r\ = £|1 — T/T c \. At T = T c = k/tt and finite 
F, Eq. ([|) is corrected by a term ~ | lnF|~^ so that mobility does not vanish ||. Also, for 
T < T c , Eq. (f|) holds only if F~\ 1 ~ T / T < : \ ^> 1. If this condition does not hold, the original 
Eq. (48) in Ref. || has to be used. According to RG analysis, these results are valid close to 
T c and in the large L, small g, and small F limit with the crossover regime characterized by 
two effective lengths: L g ~ ^/ 2 I 1 ~ T / Tc l an( j L F ^, p- 1 / 2 (g — ng 2 /2T 2 is the bare coupling 
constant). 

Numerical simulations of the model @) with the Hamiltonian @) in the limit of small 



g were recently performed by Batrouni and Hwa |L3| in the context of a randomly pinned 
planar flux array. They found no sign of phase transition in statics but in dynamics they 
observed a phase transition which is, however, only in qualitative agreement with RG pre- 
dictions (H|). The constant ( extracted from their data is significantly smaller than that in 
Eq. (U). On the other hand, the behavior of the surface under the influence of a strong 
driving force including the KPZ nonlinearity also has many interesting features and it has 
been the subject of recent investigations 

After a brief description of the simulated dynamics we will present our numerical results 
for the dynamics of the model (|I]) which maps to the model (0) with coupling constant g of 
order 1. 

Every surface configuration can be completely specified by a collection of column height 
variables C = {hi, h 2 , . . .}. The dynamics of the model is determined by the transition rates 
W(C — ► C) which specify how the system evolves from a given configuration C into a new 
configuration C The probability P(C,t) that the surface has configuration C at time t is 
determined by the following master equation in terms of these transition rates: 

a t p(c,t) = £wc" -> c)p(C,t) 

c 

-W(C -> C')P(C,t)}. (5) 

Without driving force, the system evolves to equilibrium and the transition rates satisfy 
the detailed balance condition: W(x) = W(— x)e~ x , where x = (3 AH, and AH is change of 



energy. The driving force F is included by simply adding a term AnF to AH in W, i.e., 
W = W(/3[AH + AnF]) where An is a local change in the height between configurations C 
and C. The commonly used choice of W is the Metropolis rate: 

W((3[AH + AnF}) = min{l, e -^+^nF)y ^ 

The Monte Carlo (MC) simulations presented in this work were performed on the two- 
dimensional square lattice of linear dimension L = 64, and with periodic boundary condi- 
tions. The lattice was divided into two sublattices. During the first half of the time step, all 
heights hi of one sublattice were simultaneously updated by increasing or decreasing them 
(independently) by one unit keeping the heights of the other sublattice fixed. The moves 
are then accepted or rejected according to the Metropolis rule (^) with Hamiltonian (|]) 
and constant k = 2. In the second half of the time step, the second sublattice is upgraded 
keeping the first one fixed. 

Starting with the equilibrated configurations saved after measuring the static properties 
of the system ||, the force was turned on by implementing (||). The velocity of the growing 
surface averaged over different realizations of the disorder was monitored as a function of 
MC steps. Typically, up to 10 4 initial MC steps were discarded since the system requires 
some time to reach its stationary state characterized by a uniform velocity. Measurements 
were performed over additional MC steps whose lengths depended on the values of F and 
T. The length of the runs ranged from 5 x 10 4 (for large F and large T) to 10 6 MC steps 
(for small F and small T). The average surface heights (in lattice units a) at the end 
of the MC runs were between several thousand steps for large F and T to several dozen 
steps for small F and T. For every MC run, the measurements of the surface velocity were 
grouped into several groups (usually about ten) and the average velocity of these groups 
with corresponding error bars is presented in Fig. [I] and Fig. |2|. A practical problem in 
these simulations is to measure the response of the system to very small force because very 
long MC runs are required in order to get reliable statistics. We started with the forces of 
order 1 and then decreased F gradually toward the lowest value ( F = 0.01) for which the 
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data analysis still suggests that the surface is moving with uniform velocity although the 
error bars are much larger compared with the measurements with larger F and T (see Fig. 
0) . For smaller values of F we could not extract reliable statistics within the computer time 
available. 

Figure |I| shows the behavior of the mobility, \x = v/F, as a function of temperature for 
different driving forces, F, while Fig. [| shows the log-log plot of \x versus F for different 
T. The sample averages were performed over 50 realizations of the disorder. It is clear 
from the figures that the system has a phase transition from the regime at higher T, where 
mobility is finite and independent of temperature and force, to the nonlinear regime at 
lower temperatures where \x depends on F and T. The transition itself is very broad and 
the position of the critical temperature T c is not very clear. The straight lines in Fig. ^| are 
the best fits to the fitting equations ln(v/F) = a(T) + b(T) InF. Only the six lowest values 
of F from Fig. | were included in the fit: F = 0.010, 0.015, 0.025, 0.040, 0.065, and 0.100. 
The slope of the lines, or coefficient b, corresponds to the exponent rj in formula (^). Figure 
|3] shows a comparison of the exponent rj plotted according to Eq. (TJ) (dotted line) and the 
corresponding values extracted from the data in Fig. ^| (circles). Note that at T = T c Eq. 
(||) has to be corrected with the above mentioned logarithmic contribution due to finiteness 
of F so that disagreement between the dotted line and the numerical results is expected at 
and very close to T c . Generally speaking, there is an agreement between RG calculation and 
the numerical results. 

The finite size effects are examined by repeating the simulations for a few temperatures 
below T c on the system size L = 128 and with sample averages over 25 realizations of the 
disorder. No significant difference with respect to the L = 64 results was noticed. This 
is expected since, according to the RG analysis ||, the size L = 64 is already larger than 
the effective crossover length: Lp « 1/VF = 10 (which is estimated using the smallest 
simulated value of F). 

To summarize, numerical simulations based on the Metropolis-type dynamics for which 
the local detailed balance condition is always satisfied were performed and compared with 
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predictions of RG calculations. The broad transition becomes sharper as applied force 
becomes smaller. The broadening is a consequence of the presence of the spatially varying 
pinning potential (due to disorder) with the strong coupling constant, the applied uniform 
external force, and (to a smaller degree) the finite system size. 

The comparison of the numerical results and RG predictions suggests that there is a 
qualitative and to some degree even quantitative agreement between the two, although 
numerical results suggest that the critical temperature is shifted toward higher values in 
comparison to the RG result. The deviation from linear response occurs at temperatures 
larger than T c of statics (the numerical value for T c of statics for model (0) is T c = 0.643 ± 
0.006 ||). This is probably due to the strong coupling regime where perturbative RG 
results are not expected to work. The recent self-consistent, nonperturbative Hartree type 
calculations for relaxational dynamics show that the critical temperature slowly increases 
with g if g is larger than some characteristic value below which T c does not depend on g ||15|| . 
It is believed that this discrepancy between critical temperatures in statics and dynamics is 
an effect of the existence of many metastable states, and it is also found in other physical 
systems 
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Foundation under Grant No. PHY89-04035. 



7 



REFERENCES 



[1] H. van Beijeren and I. Nolden, in Structure and Dynamics of Surfaces II, edited by W. 
Schommers and P. von Blanckenhagen, Springer- Verlag, Berlin 1987. 

[2] J. Krug and H. Spohn, in Solids far from Equilibrium, edited by C. Godreche, Cambridge 
University Press, Cambridge, England, 1991. 

[3] Kinetics of Ordering and Growth at Surfaces, edited by M. Lagally, Plenum, New York 
1990; Dynamic of Fractal Surfaces, edited by F. Family and T. Vicsek, World Scientific, 
Singapore, 1991. 

[4] S. T. Chui and J. D. Weeks, Phys. Rev. B 14, 4978 (1976); Phys. Rev. Lett. 40, 733 
(1978). 

[5] P. Nozieres and F. Gallet, J. Phys. (Paris) 48, 353 (1987); P. Nozieres, in Solids far from 
Equilibrium, edited by C. Godreche, Cambridge University Press, Cambridge, England, 
1991. 

[6] Y.-C. Tsai and Y. Shapir, Phys. Rev. Lett. 69, 1773 (1992); Phys. Rev. E 50, 3546 
(1994). 

[7] J. Toner and D. P. DiVincenzo, Phys. Rev. B 41, 632 (1990). 

[8] D. Cule and Y. Shapir, Phys. Rev. Lett. 74, 114 (1995); Phys. Rev. B (to be published). 

[9] M. P. A. Fisher, Phys. Rev. Lett. 62, 1415 (1989); T. Nattermann, I. Lyuksyutov and 
M. Schwartz, Europhys. Lett. 16, 295 (1991). 

[10] W. J. Shugard, J. D. Weeks, and G. H. Gilmer, Phys. Rev. Lett. 41, 1399 (1978). 

[11] F. Gallet, S. Balibar, and E. Rolley, J. Phys. (Paris) 48, 369 (1987). 

[12] Y. Y. Goldschmidt and B. Schaub, Nucl. Phys. B 251, 77 (1985). 

[13] G. G. Batrouni and T. Hwa, Phys. Rev. Lett. 72, 4133 (1994). 



8 



[14] J. Krug (unpublished); T. Hwa, M. Kardar, and M. Paczuski, Phys. Rev. Lett. 66 441, 
(1990); M. Rost and H. Spohn, Phys. Rev. E 49, 3709 (1994); Y.-C. Tsai and Y. Shapir, 
ibid. 50, 4445 (1994). 

[15] D. Cule and Y. Shapir, Phys. Rev. B 51, 3305 (1995). 

[16] H. Kinzelbach and H. Horner, J. Phys. (France) I 3, 1329 (1993); H. Horner, Z. Phys. 
B 86, 291 (1992); A. Crisanti, H. Horner, and H. J. Sommers, ibid. 92, 257 (1993). 



9 



FIGURES 

FIG. 1. The dependence of the mobility \x = v/F as function of temperature for different F. 
The system size is L = 64 and sample averaging is performed over 50 realizations of disorder. The 
full curves are the guides to the eye. 

FIG. 2. The log-log plot of the mobility versus driving force. The straight lines are the best 
fits to ln(v/F) = a + b\n(F) including only six the smallest values of F for each T in Fig. [l|. 

FIG. 3. Plot of the coefficient b(T) (circles) from the fitting equation (see text and Fig. |2|). 
The dotted line is the analytical prediction (0). 
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